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Preface 



Over the years we have had the rare opportunity to teach small classes of 
intelligent and strongly motivated economics students who found nonlinear 
dynamics inspiring and wanted to know more. This book began as an 
attempt to organise our own ideas on the subject and give the students 
a fairly comprehensive but reasonably short introduction to the relevant 
theory. Cambridge University Press thought that the results of our efforts 
might have a more general audience. 

The theory of nonlinear dynamical systems is technically difficult and 
includes complementary ideas and methods from many different fields of 
mathematics. Moreover, as is often the case for a relatively new and fast 
growing area of research, coordination between the different parts of the 
theory is still incomplete, in spite of several excellent monographs on the 
subject. Certain books focus on the geometrical or topological aspects of 
dynamical systems, others emphasise their ergodic or probabilistic proper- 
ties. Even a cursory perusal of some of these books will show very significant 
differences not only in the choice of content, but also in the characterisa- 
tions of some fundamental concepts. (This is notoriously the case for the 
concept of attractor.) 

For all these reasons, any introduction to this beautiful and intellectually 
challenging subject encounters substantial difficulties, especially for non- 
mathematicians, as are the authors and the intended readers of this book. 
We shall be satisfied if the book were to serve as an access to the basic con- 
cepts of nonlinear dynamics and thereby stimulate interest on the part of 
students and researchers, in the physical as well as the social sciences, with 
a basic mathematical background and a good deal of intellectual curiosity. 

The book includes those results in dynamical system theory that we 
deemed most relevant for applications, often accompanied by a common- 
sense interpretation. We have also tried, when necessary, to eliminate the 
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xii 

confusion arising from the lack of consistent and universally accepted defi- 
nitions of some concepts. Full mathematical proofs are usually omitted and 
the reader is referred either to the original sources or to some more recent 
version of the proofs (with the exception of some ‘canonical’ theorems whose 
discussion can be found virtually in any textbook on the subject). 

We devote an unusually large space to the discussion of stability, a subject 
that in the past played a central role in the theory of differential equations 
and related applied research. The fundamental monographs on stability 
were published in the 1960s or early 1970s yet there is surprisingly little 
reference to them in modern contemporary research on dynamical systems. 
We have tried to establish a connection between the classical theory of 
stability and the more recent discussions of attracting sets and attractors. 

Although the word ‘chaos’ does not appear in the title, we have dedicated 
substantial attention to chaotic sets and attractors as well as to ‘routes to 
chaos’. Moreover, the geometric or topological properties of chaotic dynam- 
ics are compared to their measure-theoretic counterparts. 

We provide precise definitions of some basic notions such as neighbour- 
hood, boundary, closure, interior, dense set and so on, which mathemati- 
cians might find superfluous but, we hope, will be appreciated by students 
from other fields. 

At an early stage in the preparation of this book, we came to the conclu- 
sion that, within the page limit agreed upon with the publisher, we could not 
cover both theory and applications. We squarely opted for the former. The 
few applications discussed in detail belong to economics where our compara- 
tive advantages lie, but we emphasised the multi-purpose techniques rather 
than the specificities of the selected models. 

The book includes about one hundred exercises, most of them easy and 
requiring only a short time to solve. 

In 1992, Cambridge University Press published a book on Chaotic Dy- 
namics by the first author, which contained the basic concepts of chaos the- 
ory necessary to perform and understand numerical simulations of difference 
and differential equations. That book included a user-friendly software pro- 
gram called DMC (Dynamical Models Cruncher). A refurbished, enlarged 
and Windows-compatible version of the program is available, at no cost, 
from the webpage 

<http://uk.cambridge.org/economics/catalogue/0521558743> 
along with a workbook of computer exercises coordinated with the ‘paper 
and pencil’ exercises found in the book. The webpage will also be used to 
circulate extra exercises, selected solutions and, we hope, comments and 
criticisms by readers. 
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Statics and dynamics: some elementary concepts 



Dynamics is the study of the movement through time of variables such as 
heartbeat, temperature, species population, voltage, production, employ- 
ment, prices and so forth. 

This is often achieved by means of equations linking the values of vari- 
ables at different, uniformly spaced instants of time, i.e., difference equa- 
tions, or by systems relating the values of variables to their time derivatives, 
i.e., ordinary differential equations. Dynamical phenomena can also be 
investigated by other types of mathematical representations, such as par- 
tial differential equations, lattice maps or cellular automata. In this book, 
however, we shall concentrate on the study of systems of difference and 
differential equations and their dynamical behaviour. 

In the following chapters we shall occasionally use models drawn from 
economics to illustrate the main concepts and methods. However, in general, 
the mathematical properties of equations will be discussed independently of 
their applications. 



1.1 A static problem 

To provide a first, broad idea of the problems posed by dynamic vis-a-vis 
static analysis, we shall now introduce an elementary model that could be 
labelled as ‘supply-demand-price interaction in a single market’. Our model 
considers the quantities supplied and demanded of a single good, defined as 
functions of a single variable, its price, p. In economic parlance, this would 
be called partial analysis since the effect of prices and quantities determined 
in the markets of all other goods is neglected. It is assumed that the demand 
function D(p) is decreasing in p (the lower the price, the greater the amount 
that people wish to buy), while the supply function S(p) is increasing in p 
(the higher the price, the greater the amount that people wish to supply). 
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Fig. 1.1 The static partial equilibrium model 



For example, in the simpler, linear case, we have: 



D(p ) = a — bp 
S(p) = —m + sp 



( 1 . 1 ) 



and a, b, m and s are positive constants. Only nonnegative values of these 
variables are economically meaningful, thus we only consider D,S,p > 0. 
The economic equilibrium condition requires that the market of the 
good clears, that is demand equals supply, namely: 



D(p) = S(p) 



( 1 . 2 ) 



or 



a — bp = — m + sp. 



static solution Mathematically, the solution to our problem is the value 
of the variable that solves (1.2) (in this particular case, a linear equation). 
Solving (1.2) for p we find: 

a + m 



where p is usually called the equilibrium price (see figure l.l ). 1 We 
call the problem static since no reference is made to time or, if you prefer, 

1 11 if demand curve D' in figure 1.1 is provided to make the point that, with no further con- 
straints on parameter values, the equilibrium price could imply negative equilibrium quantities 
of supply and demand. To eliminate this possibility we further assume that 0 < m/s < a/b , as 
is the case for the demand curve D. 




1.2 A discrete-time dynamic problem 
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everything happens at the same time. Notice that, even though the static 
model allows us to find the equilibrium price of the good, it tells us nothing 
about what happens if the actual price is different from its equilibrium value. 



1.2 A discrete-time dynamic problem 

The introduction of dynamics into the model requires that we replace the 
equilibrium condition (1.2) with some hypothesis concerning the behaviour 
of the system off-equilibrium, i.e. , when demand and supply are not equal. 
For this purpose, we assume the most obvious mechanism of price adjust- 
ment: over a certain interval of time, the price increases or decreases in 
proportion to the excess of demand over supply, ( D — S ) (for short, excess 
demand). Of course, excess demand can be a positive or a negative quan- 
tity. Unless the adjustment is assumed to be instantaneous, prices must 
now be dated and p n denotes the price of the good at time n, time being 
measured at equal intervals of length h. Formally, we have 

Pn+h = Pn + hd[D(p n ) - S(p n )}. (1.3) 

Since h is the period of time over which the adjustment takes place, 6 can 
be taken as a measure of the speed of price response to excess demand. For 
simplicity, let us choose h — 1, 0 = 1. Then we have, making use of the 
demand and supply functions (1.1), 

p n+ i = a + m+ (1 - b- s)p n . (1.4) 

In general, a solution of (1.4) is a function of time p(n ) (with n taking 
discrete, integer values) that satisfies (1.4). 2 

DYNAMIC solution To obtain the full dynamic solution of (1.4), we begin 
by setting a = a + m, (3 = (1 — b — s) to obtain 

Pn+i = a + f3p n . (1.5) 

To solve (1.5), we first set it in a canonical form, with all time-referenced 
terms of the variable on the left hand side (LHS), and all constants on the 
right hand side (RHS), thus: 

p n+ ±-(3p n = a. (1.6) 

Then we proceed in steps as follows. 

2 We use the forms p n and p(n) interchangeably, choosing the latter whenever we prefer to 
emphasise that p is a function of n. 
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STEP 1 We solve the homogeneous equation, which is formed by setting 
the RHS of (1.6) equal to 0, namely: 

p n+ i-(3p n = 0. (1.7) 

It is easy to see that a function of time p(n) satisfying (1.7) is p{n) = Cf3 n , 
with C an arbitrary constant. Indeed, substituting in (1.7), we have 

C(3 n+l - /3C/3 n = C(3 n+1 - C(3 n+l = 0. 



STEP 2 We find a particular solution of (1.6), assuming that it has a 
form similar to the RHS in the general form. Since the latter is a constant, 
set p(n) = k, k a constant, and substitute it into (1.6), obtaining 

k — /3k = a 



so that 



k 



a 



1-/9 



a + m 
b + s 



= p 



again! 



It follows that the p(ri) = p is a solution to (1.6) and the constant (or 
stationary) solution of the dynamic problem is simply the solution of the 
static problem of section 1.1. 



STEP 3 Since (1.6) is linear, the sum of the homogeneous and the particular 
solution is again a solution, 3 called the general solution. This can be 
written as 

p(n)=p+C/3 n . (1.8) 

The arbitrary constant C can now be expressed in terms of the initial con- 
dition. Putting p(0) = po, and solving (1.8) for C we have 

Po = p + C/3° =p + C 

whence C = Po~P, that is, the difference between the initial and equilibrium 
values of p. The general solution can now be re-written as 

p(n) =p+ (po ~p)P n . (1.9) 

Letting n take integer values 1,2,.. ., from (1.9) we can generate a sequence 
of values of p, a ‘history’ of that variable (and consequently, a history of 
quantities demanded and supplied at the various prices), once its value at 
any arbitrary instant of time is given. Notice that, since the function p n +i = 



3 This is called the superposition principle and is discussed in detail in chapter 2 section 2.1. 




1.2 A discrete-time dynamic problem 
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f(p n ) is invertible, i.e., the function f~ L is well defined, p n -\ = / _1 (p n ) 
also describes the past history of p. 

The value of p at each instant of time is equal to the sum of the equilib- 
rium value (the solution to the static problem which is also the particular, 
stationary solution) and the initial disequilibrium (po — p), amplified or 
dampened by a factor f3 n . There are therefore two basic cases: 

(i) \/3\ > 1. Any nonzero deviation from equilibrium is amplified in time, 
the equilibrium solution is unstable and as n — * +oo, p n asymptotically 
tends to +oo or — oo. 

(ii) \(3\ < 1. Any nonzero deviation is asymptotically reduced to zero, 
p n — > p as n — > Too and the equilibrium solution is consequently stable. 

First-order, discrete-time equations (where the order is determined as the 
difference between the extreme time indices) can also have fluctuating be- 
haviour, called improper oscillations , 4 owing to the fact that if (3 < 0, 
/ 3 n will be positive or negative according to whether n is even or odd. Thus 
the sign of the adjusting component of the solution, the second term of the 
RHS of (1.9), oscillates accordingly. Improper oscillations are dampened if 
(3 > — 1 and explosive if /3 < — 1. 

In figure 1.2 we have two representations of the motion of p through time. 
In figure 1.2(a) we have a line defined by the solution equation (1.5), and the 
bisector passing through the origin which satisfies the equation p n + i = p n - 
The intersection of the two lines corresponds to the constant, equilibrium 
solution. To describe the off-equilibrium dynamics of p, we start on the 
abscissa from an initial value po y p. To find p\ , we move vertically to 
the solution line and sidewise horizontally to the ordinate. To find p- 2 , 
we first reflect the value of p± by moving horizontally to the bisector and 
then vertically to the abscissa. From the point p\ , we repeat the procedure 
proposed for po (up to the solution line, left to the ordinate), and so on and 
so forth. The procedure can be simplified by omitting the intermediate step 
and simply moving up to the solution line and sidewise to the bisector, up 
again, and so on, as indicated in figure 1.2(a). It is obvious that for |/?| < 1, 
at each iteration of the procedure the initial deviation from equilibrium 
is diminished again, see figure 1.2(b). For example, if (3 = 0.7, we have 
/ 3 2 = 0.49, (3 3 = 0.34, . . . ,/3 10 s=a 0.03, . . .) and the equilibrium solution is 
approached asymptotically. 

The reader will notice that stability of the system and the possibility 

4 The term improper refers to the fact that in this case oscillations of variables have a ‘kinky’ 
form that does not properly describe the smoother ups and downs of real variables. We discuss 
proper oscillations in chapter 3. 
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Fig. 1.2 Convergence to p in the discrete-time partial equilibrium model 



of oscillatory behaviour depends entirely on (3 and therefore on the two 
parameters b and s, these latter denoting respectively the slopes of the 
demand and supply curves. The other two parameters of the system, a and 
m, determine a and consequently they affect only the equilibrium value 
p. We can therefore completely describe the dynamic characteristics of the 
solution (1.9) over the parameter space (6, s). The boundary between stable 
and unstable behaviour is given by \/3\ = 1, and convergence to equilibrium 
is guaranteed for 

- 1 < (3 < 1 
2 > b + s > 0. 

The assumptions on the demand and supply functions imply that 6, s > 0. 
Therefore, the stability condition is [b + s) < 2, the stability boundary is 
the line (6 + s) = 2, as represented in figure 1.3. Next, we define the curve 
f3 = 1 — (6 + s) = 0| separating the zone of monotonic behaviour from 
that of improper oscillations, which is also represented in figure 1.3. Three 
zones are labelled according to the different types of dynamic behaviour, 
namely: convergent and monotonic; convergent and oscillatory; divergent 
and oscillatory. Since b, s > 0, we never have the case (3 > 1, corresponding 
to divergent, nonoscillatory behaviour. 

If \(3\ > 1 any initial difference (po — p) is amplified at each step. In this 
model, we can have \(3\ > 1 if and only if (3 < —1. Instability, then, is due 
to overshooting. Any time the actual price is, say, too low and there is 
positive excess demand, the adjustment mechanism generates a change in 
the price in the ‘right’ direction (the price rises) but the change is too large. 



1.3 A continuous-time dynamic problem 
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Fig. 1.3 Parameter space for the discrete-time partial equilibrium model 



After the correction, the new price is too high (negative excess demand) and 
the discrepancy from equilibrium is larger than before. A second adjustment 
follows, leading to another price that is far too low, and so on. We leave 
further study of this case to the exercises at the end of this chapter. 



1.3 A continuous-time dynamic problem 

We now discuss our simple dynamical model in a continuous-time setting. 
Let us consider, again, the price adjustment equation (1.3) (with 0 = 1, 
h > 0) and let us adopt the notation p(n ) so that 

p{n + h ) = p(n) + h ( D[p(n)\ — 5[p(?r)]) . 

Dividing this equation throughout by h, we obtain 

P(n + h) — p(n) 

- = D[p(n)\ - S[p(n)] 

whence, taking the limit of the LHS as h — » 0, and recalling the definition 
of a derivative, we can write 

d ^ = D\p(n)]-S[p(n)}. 

Taking the interval h to zero is tantamount to postulating that time is a 
continuous variable. To signal that time is being modelled differently we 
substitute the time variable ngZ with i £ 1 and denote the value of p at 
time t simply by p, using the extended form p(t) when we want to emphasise 
that price is a function of time. We also make use of the efficient Newtonian 
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notation dx(t)/dt = x to write the price adjustment mechanism as 

^ = p = D(p) - S{p) = (a + m) - (6 + s)p. (1-10) 

Equation (1.10) is an ordinary differential equation relating the values of 
the variable p at a given time t to its first derivative with respect to time 
at the same moment. It is ordinary because the solution p(t) is a function 
of a single independent variable, time. Partial differential equations, whose 
solutions are functions of more than one independent variable, will not be 
treated in this book, and when we refer to differential equations we mean 
ordinary differential equations. 



DYNAMIC solution The dynamic problem is once again that of finding a 
function of time p(t) such that (1.10) is satisfied for an arbitrary initial 
condition p(0) = po. 

As in the discrete-time case, we begin by setting the equation in canonical 
form, with all terms involving the variable or its time derivatives on the LHS, 
and all constants or functions of time (if they exist) on the RHS, thus 

p + (b + s)p = a + m. (1-11) 

Then we proceed in steps as follows. 



STEP 1 We solve the homogeneous equation, formed by setting the RHS of 
(1.11) equal to 0, and obtain 

p + (b + s)p = 0 or p = —(6 + s)p. (1-12) 

If we now integrate (1.12) by separating variables, we have 

/f = - < '>+'> ) / * 

whence 

lnp(f) = — (6 + s)t + A 

where A is an arbitrary integration constant. Taking now the antilogarithm 
of both sides and setting e A = C , we obtain 

p(t) = Ce~^ b+s)t . 




1.3 A continuous-time dynamic problem 
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STEP 2 We look for a particular solution to the nonhomogeneous equation 
(1.11). The RHS is a constant so we try p = k, k a constant and conse- 
quently p = 0. Therefore, we have 

p = 0 = (a + m) — (b + s)k 



whence 



k 



a + m 
b + s 



= p. 



Once again the solution to the static problem turns out to be a special 
(stationary) solution to the corresponding dynamic problem. 



STEP 3 Since (1.12) is linear, the general solution can be found by summing 
the particular solution and the solution to the homogeneous equation, thus 

p(t) = p T Ce~ ( ' b+S ^ t . 

Solving for C in terms of the initial condition, we find 
p(0) =p 0 =p + C and C=(p 0 -p). 

Finally, the complete solution to (1.10) in terms of time, parameters, initial 
and equilibrium values is 

Pit) =P+ (po -p)e~( b+s)t . (1.13) 

As in the discrete case, the solution (1.13) can be interpreted as the sum 
of the equilibrium value and the initial deviation of the price variable from 
equilibrium, amplified or dampened by the term e -( b + s ) t . Notice that in 
the continuous-time case, a solution to a differential equation p = f(p) 
always determines both the future and the past history of the variable p, 
independently of whether the function / is invertible or not. In general, we 
can have two main cases, namely: 

(i) (6 + s) >0 Deviations from equilibrium tend asymptotically to zero as 
t — * Too. 

(ii) (6 T s) < 0 Deviations become indefinitely large as t — > Too (or, equiv- 
alently, deviations tend to zero as t — > — oo). 

Given the assumptions on the demand and supply functions, and therefore 
on b and s, the explosive case is excluded for this model. If the initial price 
is below its equilibrium value, the adjustment process ensures that the price 
increases towards it, if the initial price is above equilibrium, the price de- 
clines to it. (There can be no overshooting in the continuous-time case.) In 
a manner analogous to the procedure for difference equations, the equilibria 
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Fig. 1.4 The continuous-time partial equilibrium model 



of differential equations can be determined graphically in the plane (p,p) 
as suggested in figure 1.4(a). Equilibria are found at points of intersection 
of the line defined by (1.10) and the abscissa, where p = 0. Convergence 
to equilibrium from an initial value different from the equilibrium value is 
shown in figure 1.4(b). 

Is convergence likely for more general economic models of price adjust- 
ment, where other goods and income as well as substitution effects are taken 
into consideration? A comprehensive discussion of these and other related 
microeconomic issues is out of the question in this book. However, in the 
appendixes to chapter 3, which are devoted to a more systematic study of 
stability in economic models, we shall take up again the question of conver- 
gence to or divergence from economic equilibrium. 

We would like to emphasise once again the difference between the discrete- 
time and the continuous-time formulation of a seemingly identical problem, 
represented by the two equations 



Pn+i - p n = (a + m) - (6 + s)p n (1.4) 

p = (a + m) — (b + s)p. (1.10) 

Whereas in the latter case (6+s) > 0 is a sufficient (and necessary) condition 
for convergence to equilibrium, stability of (1.4) requires that 0 < (b+s) < 2, 
a tighter condition. 

This simple fact should make the reader aware that a naive translation of 
a model from discrete to continuous time or vice versa may have unsuspected 
consequences for the dynamical behaviour of the solutions. 
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To move from the elementary ideas and examples considered so far to a 
more general and systematic treatment of the matter, we need an appro- 
priate mathematical framework, which we introduce in this section. When 
necessary, the most important ideas and methods will be discussed in greater 
detail in the following chapters. For the sake of presentation, we shall begin 
with continuous-time systems of differential equations, which typically take 
the canonical form 



- = x = f{x) (1-14) 

where / is a function with domain U , an open subset of M m , and range 
M m (denoted by f: U — * M m ). The vector 5 x = (aq, X2, ■ ■ ■ , x m ) T denotes 
the physical variables to be studied, or some appropriate transformations 
of them; i 6 M indicates time. The variables Xi are sometimes called ‘de- 
pendent variables’ whereas t is called the ‘independent variable’. 

Equation (1.14) is called autonomous when the function / does not 
depend on time directly, but only through the state variable x. In this 
book we shall be mainly concerned with this type of equation, but in our 
discussions of stability in chapters 3 and 4 we shall have something to say 
about nonautonomous equations as well. 

The space M m , or an appropriate subspace of dependent variables — that 
is, variables whose values specify the state of the system — is referred to 
as the state space. It is also known as the phase space or, sometimes, 
the configuration space, but we will use only the first term. Although 
for most of the problems encountered in this book the state space is the 
Euclidean space, we occasionally discuss dynamical systems different from 
M m , such as the unit circle. The circle is a one-dimensional object embedded 
in a two-dimensional Euclidean space. It is perhaps the simplest example 
of a kind of set called manifold. Roughly speaking, a manifold is a set 
which locally looks like a piece of R m . A more precise definition is deferred 
to appendix C of chapter 3, p. 98. 

In simple, low-dinrensional graphical representations of the state space 
the direction of motion through time is usually indicated by arrows pointing 
to the future. The enlarged space in which the time variable is explicitly 



Recall that the transposition operator, or transpose, designated by T , when applied to a 
row vector, returns a column vector and vice versa. When applied to a matrix, the operator 
interchanges rows and columns. Unless otherwise indicated, vectors are column vectors. 




12 



Statics and dynamics: some elementary concepts 



considered is called the space of motions. Schematically, we have 



M 


x R m 


= M 1+m 


1 


1 


1 


time 


state 


space of 




space 


motions 



The function / defining the differential equation (1.14) is also called a 
vector field, because it assigns to each point x G U a velocity vector f(x). 
A solution of (1.14) is often written as a function x(t), where x : I — > M m 
and I is an interval of R. If we want to specifically emphasise the solution 
that, at the initial time to, passes through the initial point xq, we can write 
x(t\ to, xo), where x(to\ to, Xq ) = x(to) = xo . We follow the practice of setting 
to = 0 when dealing with autonomous systems whose dynamical properties 
do not depend on the choice of initial time. 



remark 1.1 In applications, we sometimes encounter differential equations 
of the form 



d m x / dx d m x x\ 
dt™ “ V ’ df ’ " ' ’ dt m_1 ) 






(1.15) 



where d k x/dt k denotes the /cth derivative of x with respect to time. Equa- 
tion (1.15) is an autonomous, ordinary differential equation of order m, 
where m is the highest order of differentiation with respect to time appear- 
ing in the equation. It can always be put into the canonical form (1.14) by 
introducing an appropriate number of auxiliary variables. Specifically, put 



d k x 

= Zk+i, 0 < k < m — 1 
(where, for k = 0, d k x/dt k = x) so that 



z k = z k+ 1 , 1 < k < m — 1 

z m — , . . • , z m ). 

If we now denote by z € R m the m-dimensional vector (z i, . . . , z m ) T we can 
write: 

z = f(z) 

where f(z) = [z 2 , ■ ■ ■ , z m , F(z 1 , . . . , z m )] T . 



We can also think of solutions of differential equations in a different man- 
ner which is now prevalent in dynamical system theory and will be very 
helpful for understanding some of the concepts discussed in the following 
chapters. 
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If we denote by 4>t(x) = (j>(t,x) the state in M m reached by the system 
at time t starting from x , then the totality of solutions of (1.14) can be 
represented by a one-parameter family of maps 6 <j>t : U K m satisfying 



S Wt ' x)l 



t=T 



f[<f>(T,x)\ 



for all x £ U and for all r 6 I for which the solution is defined. 

The family of maps <pt{x) = <p(t,x ) is called the flow (or the flow map) 
generated by the vector field /. If / is continuously differentiable (that 
is, if all the functions in the vector are continuously differentiable), then 
for any point x'o in the domain U there exists a S(x o) > 0 such that the 
solution cj)(t, xq) through that point exists and is unique for \t\ < 6 . The 
existence and uniqueness result is local in time in the sense that 5 need not 
extend to (plus or minus) infinity and certain vector fields have solutions 
that ‘explode’ in finite time (see exercise 1.8(c) at the end of the chapter). 

When the solution of a system of differential equations x = f(x) is not 
defined for all time, a new system x = g(x) can be determined which has 
the same forward and backward orbits in the state space and such that each 
orbit is defined for all time. If ip(t, x) is the flow generated by the vector field 
g, the relation between ip and the flow (p generated by / is the following: 



ip(t,x) = <p[r(t,x),x] x G U 



and 



is a time-reparametrisation function monotonically increasing in t for all 
x £ U. 



EXAMPLE Suppose we have a system 



X = f(x) (1.16) 

with / : R m — > R m , a continuously differentiable function with flow 4>{t,x) 
defined on a maximal time interval —oo<a<0<b< +oo. Then the 

6 The terms map or mapping indicate a function. In this case, we speak of y = f(x ) as the 
image of x under the map /. If / is invertible, we can define the inverse function /“ 1 , that 
is, the function satisfying / -1 [/(#)] = x for all x in the domain of / and f[f~ 1 {y)\ = y for all 
y in the domain of /“ 1 . Even if / is not invertible, the notation f~ 1 (y) makes sense: it is the 
set of pre-images of y, that is, all points x such that f(x) = y. The terms map, mapping 
are especially common in the theory of dynamical systems where iterates of a map are used to 
describe the evolution of a variable in time. 
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(a) (b) 



Fig. 1.5 A damped oscillator in M 2 : (a) space of motions; (b) state space 

differential equation x = g{x) with g(x) : R m — * M m and 

_ /( x ) 
a{) i + ll/MII 

(where |j • |j denotes Euclidean norm), defines a dynamical system whose 
forward and backward orbits are the same as those of (1.16) but whose 
solutions are defined for all time.' 

The set of points {4>(t, xo) | t £ 1} defines an orbit of (1.14), starting from 
a given point xq. It is a solution curve in the state space, parametrised by 
time. The set {[i, 0(t,x o)] | t G 1} is a trajectory of (1.14) and it evolves 
in the space of motions. However, in applications, the terms orbit and 
trajectory are often used as synonyms. A simple example of a trajectory in 
the space of motions R x R 2 and the corresponding orbit in the state space 
M 2 is given in figure 1.5. Clearly the orbit is obtained by projecting the 
trajectory onto the state space. 

The flows generated by vector fields form a very important subset of a 
more general class of maps, characterised by the following definition. 

definition 1.1 A flow is a map cf : I ClxX->I where X is a metric 
space, that is, a space endowed with a distance function, and <j> has the 
following properties 

(a) 0(0, x) = x for every x € X (identity axiom); 

(b) <f)(t + s,x) = (j>[s, (f[t., x)] = 4>[t,(j)(s,x)\ = (j>(s + t,x), that is, time- 
translated solutions remain solutions; 

'For details see Bhatia and Szego (1970), p. 78; Robinson (1999), pp. 146-7. 
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(c) for fixed t, fit is a homeomorphism on X. 

Alternatively, and equivalently, a flow may be defined as a one-parameter 
family of maps fit : X —> X such that the properties (a)-(c) above hold for 
all t,s€R. 

remark 1.2 A distance on a space X (or, a metric on X) is a function 
X x X satisfying the following properties for all x,y G X: 

(1) d(x, y) > 0 and d(x, y) = 0 if and only if x = y. 

(2) d(x,y) = d(y,x) (symmetry); 

(3) d(x, y) < d(x, z ) + d(z, y) (triangle inequality). 

Notice that there also exist notions of distance which are perfectly mean- 
ingful but do not satisfy the definition above and therefore do not define a 
metric, for example: 

the distance between a point and a set A; 

d(x,A) = inf d(x,y). 

y&A 

the distance between two sets A and B 

d(A,B) = inf inf d(x,y). 

x^lA y(zB 

Neither of these cases satisfies property (1) in remark 1.2. However, there 
exists a ‘true’ distance between two sets which is a metric in the space of 
nonempty, compact sets, i.e., the Hausdorff distance. 8 

In this book we are mainly concerned with applications for which fi is a 
flow generated by a system of differential equations and the state space is 
an Euclidean space or, sometimes, a manifold. However, some concepts and 
results in later chapters of the book will be formulated more generally in 
terms of flows on metric spaces. 

Consider now a system of nonautonomous differential equations such that 

x = f(t,x ) (1.17) 

where f : M. x U —* M m , and assume that a unique solution exists for all 
(to,xo) €E R x U. Then we can represent solutions of (1.17) by means of 
a flow fi : M x X — > X, where X C (R x M m ). This suggests that a 



See, for example, Edgar (1990), pp. 65-6. 
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nonautonomous system x = f(t, x ) can be transformed into an equivalent 
autonomous system by introducing an arbitrary variable 9 = t and writing 

^ ~~ 1 (1.18) 
x = f(9, x). 

Notice that, by definition, the extended autonomous system (1.18) has no 
equilibrium point in X. However, if the original, nonautonomous system 
(1.17) has a uniformly stable ( uniformly , asymptotically stable) equilibrium 
point, then for the extended autonomous system (1.17), the t-axis is a stable 
(asymptotically stable) invariant set. The precise meaning of (asymptotic, 
uniform) stability will be discussed in chapters 3 and 4. 

Solutions of system (1.14) can be written in either the simpler form x(t), 
x : I — * K m , or <f>t(x) : U — > M m , or again (f(t,x), 4> : I x U — > R m , 

depending on what aspect of solutions one wants to emphasise. The notation 

4>t(x) is especially suitable for discussing discrete-time maps derived from 
continuous-time systems. 

If time t is allowed to take only uniformly sampled, discrete values, sep- 
arated by a fixed interval r, from a continuous-time flow we can derive a 
discrete-time map (a difference equation) 

Xn+r — G{,Xn) (1.19) 

where G = <fi T . Certain properties of continuous-time dynamical systems 
are preserved by this transformation and can be studied by considering the 
discrete-time systems derived from them. If the unit of measure of time is 
chosen so that r = 1, we have the canonical form 

x n+ i = G(x n ). (1-20) 

Let the symbol o denote the composition of functions, so that, / o g{x) 
means f[g(x)]. Then we write 

x n = G(x n -i) = G O G(x n - 2 ) = ... = GoGo...o G(x 0 ) = G n (x 0 ) 

where G n is the composition of G with itself n times, or the nth iteration 
of G , with n E Z + . If G is invertible and G -1 is a well defined function, G n 
with n E TL~ denotes the nth iterate of G~ x . (Note that G n (x) is not the 
nth power of G(x).) Thus, iterates of the map G (or G~ 1 ) can be used to 
determine the value of the variable x at time n, when the initial condition 
xq is fixed. 9 



9 For autonomous difference equations whose solutions do not depend on the choice of the initial 
time, in a manner analogous to our practice for autonomous differential equations, we take the 
initial time as zero. 
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remark 1.3 There exists another way of deriving a discrete-time map from 
a continuous-time dynamical system, called Poincare map, which de- 
scribes the sequence of positions of a system generated by the intersections 
of an orbit in continuous time and a given space with a lower dimension, 
called surface of section. Clearly, in this case the time intervals between 
different pairs of states of the systems need not be equal. Poincare maps 
are a powerful method of investigation of dynamical systems and we shall 
make some use of them in chapter 4, when we discuss periodic solutions and 
in chapters 7 and 8. 

Of course, there exist problems that are conceived from the beginning 
as discrete dynamical systems (difference equations). In fact, there are 
difference equations that cannot be derived from differential equations. In 
particular, this is true of noninvertible maps which have been extensively 
used in recent years in the study of dynamical problems in many applied 
fields. Intuitively, the reason why a noninvertible map cannot be a flow map 
(derived from a differential equation as explained above) is that such a map 
uniquely determines the dynamics in one time direction only whereas, under 
standard assumptions, solutions of a differential equation always determine 
the dynamics in both directions uniquely. 

remark 1.4 Orbits of differential equations are continuous curves, while 
orbits of maps are discrete sets of points. This has a number of important 
consequences, the most important of which can be appreciated intuitively. 
If the solution of an autonomous system of differential equations through a 
point is unique, two solution curves cannot intersect one another in the state 
space. It follows that, for continuous-time dynamical systems of dimension 
one and two, the orbit structure must be drastically constrained. In the 
former, simpler case, we can only have fixed points and orbits leading to (or 
away from) them; in the two-dimensional case, nothing more complex than 
periodic orbits can occur. For maps the situation is different. It remains true 
that the orbit starting from a given point in space is uniquely determined in 
the direction defined by the map. However, since discrete-time orbits, so to 
speak, can ‘jump around’, even simple, one-dimensional nonlinear maps can 
generate very complicated orbits, as we shall see in the following chapters. 

Generalising the simple examples discussed in sections 1.2 and 1.3 above, 
the stationary, equilibrium solutions of multi-dimensional dynamical sys- 
tems in both continuous and discrete time can be identified by solving sys- 
tems of equations. 

In the former case, setting x = 0 in (1.14) the set of equilibrium or 
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fixed points is defined by 

E = {x\f(x) = 0} 

that is, the set of values of x such that its rate of change in time is nil. 

Analogously, in the discrete-time case, 

Xn+l — G(x n ) 

we have 

E = {x\x — G(x ) = 0} 

that is, the set of values of x that are mapped to themselves by G. Because 
the functions / and G are generally nonlinear, there are no ready-made 
procedures to find the equilibrium solutions exactly, although geometrical 
and numerical techniques often give us all the qualitative information we 
need. Notice that linear systems typically have a unique solution, whereas 
nonlinear systems typically have either no solutions, or a finite number of 
them. It follows that only nonlinear systems may describe the interesting 
phenomenon of (finite) multiple equilibria. 

For a system of autonomous, differential equations like (1.14), a general 
solution <j>(t,x ) can seldom be written in a closed form, i.e. , as a combi- 
nation of known elementary functions (powers, exponentials, logarithms, 
sines, cosines, etc.). Unfortunately, closed-form solutions are available only 
for special cases, namely: systems of linear differential equations; one- 
dimensional differential equations (i.e., those for which m = 1); certain 
rather special classes of nonlinear differential equations of order greater than 
one (or systems of equations with m > 1). The generality of nonlinear sys- 
tems which are studied in applications escapes full analytical investigation, 
that is to say, an exact mathematical description of solution orbits cannot 
be found. Analogous difficulties arise when dynamical systems are repre- 
sented by means of nonlinear maps. In this case, too, closed-form solutions 
are generally available only for linear systems. 

The importance of this point should not be exaggerated. On the one 
hand, even when a closed-form solution exists, it may not be very use- 
ful. A handbook of mathematical formulae will typically have a hundred 
pages of integrals for specific functions, so that a given nonlinear model 
may indeed have a solution. However, that solution may not provide much 
intuition, nor much information if the solution function is not a common, 
well known function. On the other hand, in many practical cases we are 
not especially interested in determining (or approximating) exact individual 
solutions, but we want to establish the qualitative properties of an ensemble 
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of orbits starting from certain practically relevant sets of initial conditions. 
These properties can often be investigated effectively by a combination of 
mathematical, geometrical, statistical and numerical methods. Much of 
what follows is dedicated precisely to the study of some of those methods. 

Before turning to this goal, however, we review in chapter 2 the class of 
dynamical systems which is best understood: linear systems. Dynamical lin- 
ear systems in both continuous and discrete time are not terribly interesting 
per se because their behaviour is morphologically rather limited and they 
cannot be used effectively to represent cyclical or complex dynamics. How- 
ever, linear theory is an extremely useful tool in the analysis of nonlinear 
systems. For example, it can be employed to investigate qualitatively their 
local behaviour, e.g., their behaviour in a neighbourhood of a single point 
or of a periodic orbit. This is particularly important in stability analysis 
(chapters 3 and 4) and in the study of (local) bifurcations (chapter 5). 



Exercises 

1.1 Consider the discrete-time partial equilibrium model summarised in 
(1.6) given the parameter values a = 10, b = 0.2, m = 2, s = 0.1. 
Write the general solution given the initial values po = 20 and po = 
100. Calculate the values for the price at time periods 0, 1, 2, 4, 
10, 100 starting from each of the above initial values and sketch the 
trajectories for time periods 0-10. 

1.2 State a parameter configuration for the discrete-time partial equi- 
librium model that implies (3 < 0. Describe the dynamics implied 
by that choice. Using these parameter values and a = 10, m = 2, 
sketch the dynamics in the space (p n ,p n + 1 ). Draw the bisector line 
and from the chosen initial condition, iterate 3 or 4 times. Show the 
direction of movement with arrows. 

1.3 If we define the parameters as in exercise 1.1 ( b = 0.2, s = 0.1, 
a = 10, m = 2), the continuous-time, partial equilibrium model of 
(1.11) gives the constant exponent of the solution as b + s = 0.3. Let 
this be case 1. If s = 0.6, b+s = 0.8. Let this be case 2. Calculate the 
solution values for case 1 and case 2 at periods t = 0, 1, 2, 4.67, 10, 100 
starting from the initial condition po = 20. Comment on the speed 
of the adjustment process. Note the different integer values of t for 
which equilibrium in Case 2 is approximated using a precision of 1 
decimal point, 2 decimal points, 3 decimal points. 

1.4 Suppose that the good under consideration is a ‘Giffen’ good (for 
which dD/cLp > 0 and therefore b < 0). It is unlikely, but possible 
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that b+s < 0. Sketch the differential equation (1.10) under that hy- 
pothesis in the (p,p) plane, note the equilibrium point and comment 
on the adjustment process. 

Convert these higher-order differential equations to systems of first- 
order differential equations and write the resulting systems in matrix 
form: 

(a) x + x = 1 

(b) gf+0.4x-2x = 0 

(c) + 4x — 0.5x — x = 11. 

Convert the following higher-order system of differential equations 
into a system of first-order differential equations 

x + x = 1 

y-y-y = -i- 



Higher-order difference equations and systems can also be converted 
to first-order systems using auxiliary variables. A fcth-order equation 
x n+ k = G(x n+ k~ i, • • • , x n ) can be converted by setting 



Xn 



= Z 



( 1 ) 



(1) 

4+1 = X n+ l = z 



( 2 ) 



(2) 

4+1 = X n+ 2 = z 



( 3 ) 



z nll = x n+k = G(x n+k - 1, . . . , X n ) = G(4 fc) , • • • , 4 1} )- 

Convert the following difference equations into systems of first-order 
difference equations and write them in matrix form 

(a) x n+2 - ax n+ i + bx n = 1 

(b) 0.5x72+3 4“ 2x77+1 O.IX72 — 2. 

1.8 Use integration techniques to find exact solutions to the following 
differential equations and sketch trajectories where possible, assum- 
ing an initial value of x(0) = 1 

(a) x = 2 x 

( b ) * = 4 

(c) X = x 2 . 

1.9 Use the technique described in the example in section 1.4 to find 
a function g, defined over all time and such that x = g(x) has the 
same backward and forward orbits in the state space as x = x 2 . 
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1.10 Write the exact solution of the following differential equation ( Hint 
rewrite the equation as dx/dt = /xx{\ — x) and integrate, separating 
variables) and discuss the dynamics of x 

x = nx{ 1 — x) x £ [0, 1]. 




2 

Review of linear systems 



The theory of linear dynamical systems is a very well-developed area of 
research and even an introductory presentation would easily fill a large book 
in itself. Many excellent mathematical treatments of the matter exist and 
we refer the reader to them for more details. Our objective in what follows 
is to provide a review of the theory, concentrating on those aspects which 
are necessary to understand our discussion of nonlinear dynamical systems. 
We consider the general multi-dimensional case first and then, to fix certain 
basic ideas, we discuss in detail simple mathematical examples in the plane, 
as well as an extended version of the economic partial equilibrium model. 



2.1 Introduction 

Broadly speaking, we say that a phenomenon represented by a stimulus- 
response mechanism is linear if, to a given change in the intensity of the 
stimulus, there corresponds a proportional change in the response. Thus, 
postulating that saving is a linear function of income implies that a doubling 
of the level of income doubles the amount of saving as well. 

As concerns dynamical systems, we say that systems such as 



dx ■ ft \ 

Tt =x = l(x) 


x € M m iel 


(2.1) 


^n+1 = G(^X n ^) 


x 6 M m n6Z 


(2.2) 



are linear if the vector- valued functions f(x) and G(x n ) are linear according 
to the following definition: 

definition 2.1 A function f : M m — > M m is linear if f(av+/3w) = af(v) + 
/3f(w), for any a, f3 E M and v,w G W n . 
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In chapter 1 we provided a preliminary discussion of some simple dynam- 
ical systems and said that, since they were linear, the sum of solutions 
was also a solution. We can now generalise that idea and make it more 
precise, beginning with systems in continuous time. Linear systems of dif- 
ferential equations satisfy the superposition principle, defined as follows: 
if x) and f 2 (t, x ) are any two linearly independent 1 solutions of system 
(2.1) then 

S(t, x) = a(j>i(t, x) + fi(f> 2 (t, x) (2.3) 

is also a solution for any a, (3 E M. It can be verified that the superposition 
principle is valid if, and only if, the (vector- valued) function / of (2.1) is 
linear. If <f>i (t, x) and <f) 2 (t,x) are solutions to (2.1) then, respectively: 

= f[M r , x )\ 

^[< fo{t,X )] = f[<h(T,x)\. 

On the other hand, in order for (2.3) to be a solution, we must have 

( ^[‘S'(i,z)] ^ =f[S(r,x)] 
or 

a/[</>i(r, x)] + Pf[ 4 > 2 (r, x)] = /[a0i(r, x) + /3</> 2 (t, x)] 

which holds if, and only if, / is linear. 

An entirely analogous argument can be developed for the discrete-time 
system (2.2) and we leave it to the reader as an exercise. 

When / and G are vectors of linear functions, their elements can be 
written as products of constants and the variables, namely, 

fi(x) = fax 1 + f i2 x 2 H b fimXm (* = 1, 2, . . . ,.m) 

and 

Gii^x) — GilXi Gi 2 X 2 * * * ~\~ GifyiXrn (z — 1,2,..., ?7z) , 

respectively, where fi and Gi denote the ith element of / and G (called 
coordinate function), respectively, and fij,Gij (i. j = 1,2, ...,m) are 
constants. 



say that (j>\(t,x ) and 02 {t,x) are linearly independent if, and only if, a<f>i(t,x) + 
(3(j>2(t,x) = 0 for all x and all t implies that a = (3 = 0. 
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Then systems (2.1) and (2.2) can be put in a compact, matrix form, 



respectively, as 










x = Ax 


x £ R m 


(2.4) 


and 










^n+l = Bx n 


x £ M m 


(2.5) 



where A and B are m x m matrices of constants, also called the coefficient 
matrices, with typical elements fij and G t] , respectively. 

remark 2.1 Consider the equation 



x = f(x) = a + bx (2.6) 

where a and b are scalar constants. If we try definition 2.1 on the function 
f(x), letting x take on the values v and w, we have 

f(av + (3w) = a + bav + bf3w. 

But 

af(v) + (3f(w) = aa + abv + /3a + (3bw 

and, for o ^ 0, these are equal if, and only if, a + /3 = 1. Therefore the 
requirements of definition 2.1 are not satisfied and the superposition prin- 
ciple does not hold for (2.6). Functions like / are linear plus a translation 
and are called affine. 

Dynamical systems characterised by affine functions though, strictly 
speaking, nonlinear, can easily be transformed into linear systems by trans- 
lating the system without affecting its qualitative dynamics. 

We show this for the general, continuous-time case. Consider the system 

y = Ay + c (2.7) 

where y £ M m , A is an m x m matrix and c £ M m is a vector of constants. 
We can set 

x = y + k 

so that 

x = y = Ay + c = Ax — Ak + c. 

If A is nonsingular (the determinant, det(yl), is nonzero and therefore 
there are no zero eigenvalues), A -1 is well defined and we can set k = A ^c 
whence x = Ax which is linear and has the same dynamics as (2.7). Notice 
that the assumption that A is nonsingular is the same condition required 
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for determining the fixed (equilibrium) points of (2.7). Recall that a point 
y is a fixed or equilibrium point for (2.7) if y = Ay + c = 0, whence 

fj = -A^c. 

If A is singular, that is, if A has a zero eigenvalue, the inverse matrix A -1 
and, therefore, the equilibrium are indeterminate. 

The discrete-time case is entirely analogous. Consider the system 



!J a ■ l — By n T c (2.8) 

where y n E M m , B is a m x m matrix and c E M m is a vector of constants. 
We set 

Xn — Vn T k 



so that 



x n + 1 = y n + 1 + k — By n + c+ k = Bx n - Bk + c+ k. 

If we set k = — [I — B]~ 1 c (where I denotes the identity matrix), we have 
x n+ \ = Bx n , a linear system of difference equations with the same dynamics 
as (2.8). Once again k is defined if the matrix [I — B } is nonsingular, i.e. , 
if none of the eigenvalues of the matrix B is equal to one. This is also the 
condition for determinateness of the equilibrium y = [I ~ B]~ 1 c of system 
( 2 . 8 ). 

2.2 General solutions of linear systems in continuous time 

In this section we discuss the form of the solutions to the general system of 
linear differential equations (2.4) 



x = Ax x E M m . 

First, notice that the function x(t) = 0 solves (2.4) trivially. This special 
solution is called the equilibrium solution because if x = 0, x = 0 as well. 
That is, a system starting at equilibrium stays there forever. Notice that if 
A is nonsingular, x = 0 is the only equilibrium for linear systems like (2.4). 
Unless we indicate otherwise, we assume that this is the case. 

Let us now try to determine the nontrivial solutions. In chapter 1, we 
learned that, in the simple one-dimensional case, the solution has the ex- 
ponential form ke^, k and A being two scalar constants. It is natural to 
wonder whether this result can be generalised, that is, whether there exist 
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(real or complex) constant vectors u and (real or complex) constants A such 
that 

x(t) = e xt u (2.9) 

is a solution of (2.4). Differentiating (2.9) with respect to time, and substi- 
tuting into (2.4), we obtain 

A e xt u = Ae xt u 

which, for e xt ^ 0, implies 

Au = Xu or ( A — A I)u = 0 (2.10) 

where 0 is an m-dimensional null vector. A nonzero vector u satisfying 
(2.10) is called an eigenvector of matrix A associated with the eigenvalue 
A. Equation (2.10) has a nontrivial solution u ^ 0 if and only if 

det(A — XI) = 0. (2.11) 

Equation (2.11) is called the characteristic equation of system (2.4) and 
can be expressed in terms of a polynomial in A, thus 

det(A - XI) = V{X) = X m + fciA m_1 + • • • + k m - 1 A + k m = 0 (2.12) 

and V( A) is called the characteristic polynomial. The answer to our 
question is, therefore, yes, (2.9) is indeed a solution of (2.4) if A is an 
eigenvalue and u is the associated eigenvector of A. Thus, we have reduced 
the problem of solving a functional equation (2.4) to the algebraic problem 
of finding the roots of the polynomial (2.12) and solving the system of 
equations (2.10) which, for given A, is linear. Notice that for each A, if u is 
an eigenvector, so is cu for c / 0. 

remark 2.2 The coefficients ki of (2.12) depend on the elements of the 
matrix A and, in principle, can be determined in terms of those elements 
by means of rather complicated formulae. However, the first and the last, 
namely k\ and k m , can be related in a particularly simple manner, to the 
matrix A and its eigenvalues as follows: 

m 

k\ = — tr(A) = X i 

2=1 

m 

k m (— l) m det(A) = (— 1) TO n Xi 

2=1 

where tr(A) denotes the trace of A, that is, the sum of the elements on the 
main diagonal and, again, det(A) denotes the determinant of A. 
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Equation (2.12) has m roots which may be real or complex and some of 
them may be repeated. In order to render the discussion more fluid we give 
general results under the assumption that eigenvalues are distinct, followed 
by some comments on the more complicated case of repeated eigenvalues. 
As regards complex roots, if the elements of matrix A are real, complex 
eigenvalues (and eigenvectors) always appear in conjugate pairs. 

Suppose there exist n (0 < n < m) real, distinct eigenvalues A i of the 
matrix A with corresponding n real, linearly independent eigenvectors Ui, 
and p = ( m — n)/2 pairs of distinct complex conjugate 2 eigenvalues (A j, A j) 
with corresponding eigenvectors ( Uj,Uj ) where (0 < p < m/2). Then, in 
the real case, there are n linearly independent solutions of (2.4) with the 
form 

Xi(t) = e Xit Ui (i — 1,2,..., n). (2.13) 

In the complex case, each pair of eigenvalues and associated pair of eigen- 
vectors can be written, respectively, as 

(A j, A/'+i) = (-\?> -\?) = ( a j + 0j: a j ~ iftj) 

where a, (3 are real scalars, and 

(v,j,Uj+ 1 ) = ( Uj , Uj) = (aj + ibj, aj — ibj ) 

where aj and bj are m-dimensional vectors and, i 2 = — 1. The two corre- 
sponding real solutions of (2.4) are the vectors 

Xj(t) = ^(e^Uj + e^Uj) 

= e° jt aj cos ((3jt) — bj sin (/3jt) 

Xj+i(t) = - l -(e x ^Uj - 

= e ajt aj sin (f3jt) + bj cos (/3jt) 

where we have employed Euler’s formula e± l ° = cos 0 =t i sin 6 to write 

e {a±iP)t = e at [ cos(/?t ) ± isin^t)] . 

Equations (2.14) can be further simplified by transforming the real and 
imaginary parts of each element of the complex vectors ( Uj,v,j ) into the 

2 If z = a + i0 is an arbitrary complex number, then the complex number z = a — if 3 is the 
complex conjugate of z. The same notation is used for vectors. 



(2.14) 
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polar coordinates 



df* = <f cos(^°) 

bf* = Cf 1 sin(0^) l = 1, 2, . . . , m 



where ctp and are the /th elements of a,j and bj, respectively. 
Making use of the trigonometric identities 

sin(x ± y) = sin x cos y ± cos x sin y 
cos(x ± y) = cos x cos y =F sin x sin y 



(2.15) 



the solutions (2.14) of (2.4) can be written as m-dimensional vectors whose 
Ith. elements have the form 



Xj ) = Cj e ajt cos (Pjt + 4>j ) 
xf +l {t) = Cf e** 1 sm(Pjt + cj)f) 



(2.16) 



Applying the superposition principle, we can now write the general solu- 
tion of (2.4) as a function of time and m arbitrary constants, namely: 

x{t) = cixi(t) + C 2 X 2 (t) H h c m x m (t) (2.17) 



where Xi(t) (i = 1 ,...,m) are the individual solutions defined in (2.13), 
(2.14) (or (2.16)) and the m constants Ci are determined by the initial 
conditions. It is possible to show that formula (2.17) gives all of the possible 
solutions to (2.4). 

Simple inspection of (2.13) and (2.14) suggests that if any of the eigen- 
values has a positive real part, solutions initiating from a generic point (one 
for which Cj ^ 0 Vj) diverge, i.e. , variables tend to + or — oo as time tends 
to plus infinity. On the contrary, if all eigenvalues have negative real parts, 
solutions asymptotically tend to zero, that is, to the unique equilibrium 
point. Anticipating our systematic discussion of stability in chapter 3, we 
shall say that in the latter case, the equilibrium is asymptotically sta- 
ble, in the former case the equilibrium is unstable. The situation is more 
complicated if there are no eigenvalues with positive real parts but some 
eigenvalues have zero real parts. If there is only one such eigenvalue, the 
system has an attenuated form of stability — it does not converge to equi- 
librium, but it does not diverge either. (If one of the eigenvalues is zero, 
the equilibrium need not be unique.) 

Consider now, the case of a real eigenvalue Xk and set the initial conditions 
so that Ck 7 ^ 0, Ci = 0 V^. From (2.13) and (2.17) we have 

x(0) = CfcXfc(O) = c k u k ■ 



(2.18) 
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Thus, we have chosen initial conditions so as to position the system, at time 
zero, on the eigenvector associated with X k . Using (2.13), (2.17) and (2.18) 
we deduce that the solution for |f| > 0 in this case is given by 

x(t) = c k x k (t) = e Xkt c k u k = e Xkt x(0). 



Then, if the initial conditions are on the eigenvector, the system either ap- 
proaches or moves away from the equilibrium point (depending on the sign 
of the eigenvalue), along that vector. In other words, each real eigenvector 
u k defines a direction of motion in the state space that is preserved by the 
matrix A. In this case, we say that eigenvectors are invariant sets. Broadly 
speaking, a certain region of the state space is said to be invariant with 
respect to the action of a continuous- or discrete-time dynamical system, if 
an orbit starting in the set remains there forever (unless disturbed by ex- 
ogenous shocks). The speed of the motion along the eigenvector u k is given 
by x{t) = \ k e Xkt x(0) = A k x(t), that is, the speed is a constant proportion 
of the position of x. 

The complex case is more difficult to analyse, but it is possible to show 
that the plane S spanned (or generated) by the two linearly independent real 
vectors associated with a pair of complex conjugate eigenvectors ( Uj,v,j ), 



2 (uj + uj) = <ti 
-\( u 3 - %) = b j 



(2.19) 



is invariant. (Vectors aj and bj are the real and imaginary parts of (uj , Uj ) . ) 
If we choose the initial conditions on S so that the coefficients c* of the 
general solution are all zero except the two corresponding to the pair of 
complex eigenvalues (A j, A j), the orbits of the system remain on the plane 
S forever. 

Again, assuming that the matrix A has m distinct eigenvalues, we divide 
the eigenvectors (or, in the complex case, the vectors equal to the real and 
imaginary parts of them) into three groups, according to whether the cor- 
responding eigenvalues have negative, positive or zero real parts. Then the 
subsets of the state space spanned (or generated) by each group of vectors 
are known as the stable, unstable and centre eigenspaces, respectively, 
and denoted by E s , E u and E c . If m s is the number of eigenvalues with 
negative real parts, m u the number with positive real parts and m c the num- 
ber with zero real parts, the dimension of the stable, unstable and centre 
eigenspaces are m s , m u and m c , respectively, and m = m s + m u + m c . The 
subspaces E s , E u and E c are invariant. The notion of stable, unstable and 
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centre eigenspaces can be generalised to the case of repeated eigenvalues, 
but we do not discuss the question here. 

There exists a nice, compact way of writing the general solution of (2.4). 
Let X(t) be the matrix whose m columns are the m solutions defined by 
(2.13) and (2.14), or (2.16), that is, 

X(t) = [xi(t) x 2 {t) . . . x m (t)\ . 

X(t) is called the fundamental matrix solution of x = Ax and it has 
some interesting properties, in particular: 

(i) X(t) = e tA X(0), where e tA is a converging, infinite sum of (m x m) 
matrices defined as 

e tA = ^I + tA+ t -A 2 + --- + t -A n + --y, 

(ii) X(t) = AX(t), that is, X(t) satisfies (2.4) as a matrix equation; 

(iii) if the solutions Xi(t) are linearly independent, X(t) is nonsingular for 
all t. 

remark 2.3 The reader can verify, looking at (2.13) and (2.14), that the 
columns of X(0) are of two types: 

(1) for a solution of type (2.13) (real eigenvalues and eigenvectors), the 
columns are equal to m ; 

(2) for each pair of solutions of type (2.14) with (uj, Uj + \) = ( aj + ibj , aj — 
ibj ) (complex conjugate eigenvalues and eigenvectors), the correspond- 
ing two columns are equal to the real vectors aj and bj , respectively. 

From (2.17) we have 

x{t) = X(t)c 

where c = (ci, c 2 , ■ . ■ , c m ) T and therefore, from property (i) of matrix X , 
the solutions of (2.4) can be written as 

x{t) = e tA x{ 0). (2.20) 

We can verify that (2.20) is actually a solution of (2.4) by differentiating 
with respect to time and obtaining 

x(t) = Ae tA x(0) = Ax(t) 

which is (2.4). 

When there are m distinct, real eigenvalues (and m linearly independent 
real eigenvectors), the fundamental matrix solution X(t) can be used to 
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transform the matrix A into a diagonal matrix having the same eigenvalues 
as A. The reader can verify that 



AX(t) = X(t) A A 



/A, 

V 0 



Because X(t) is nonsingular we have 



X(t)~ l AX(t) = A. 




When there are distinct, complex eigenvalues the transformation, as usual, 
is more complicated and the resulting operation is a block diagonalisation. 
Consider a simple case in which A is a (2 x 2) matrix with two complex 
conjugate eigenvalues A = a + i/3, A = a — i/3, and two corresponding eigen- 
vectors u = a + ib and u = a — ib. Recalling that tr(A) = A + A = 2a and 
det(A) = AA = a 2 + (3 2 , the characteristic polynomial in this case can be 
written as 

A 2 - 2aA + (a 2 + (I 2 ) = 0. 

Making use of a theorem by Cayley and Hamilton stating that if V(X) is 
the characteristic polynomial of a matrix A, then the matrix polynomial 
V(A) = 0 (where 0 is a matrix of the same order as A in which all elements 
are null), we can write 

A 2 - 2 aA + (a 2 + j3 2 )! = 0 



or, factoring, 

(A - XI) (A - XI) = (A - a/) 2 + /3 2 I = 0. (2.21) 

If x is any nonzero real (n x 1) vector and y = [A — al)x/ (3 ( x , y are linearly 
independent if fl yf 0), we have 

Ax = (3y + ax. (2.22) 

Pre-multiplying (2.22) by (A — al) and making use of (2.21) we have 



Then defining 



Ay = ay — (3x. 



P = 



x (!) y { !) 

x ( 2 ) y ( 2 ) 



where x^\ y® denote the ith elements of the vectors x,y, respectively, we 
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have 



or 



AP = P 



a 

P 




P~ l AP = 



a —f3 
(3 a 



(2.23) 



The matrix on the RHS of (2.23) is obviously similar to A (they have the 
same eigenvalues). In the multi-dimensional case, with real and complex 
conjugate eigenvalues and eigenvectors (but no multiple eigenvalues), the 
matrix of coefficients A can be transformed into a similar block diagonal 
matrix with the general form 



... o\ 

A* 

Bj 

Vo ••• 



where A j are real eigenvalues and Bj are (2 x 2) blocks of the same form as 
(2.23). 

We conclude this section with a brief comment on the complications aris- 
ing from repeated or multiple eigenvalues. Consider a simple system of two 
differential equations for which there is a double real eigenvalue Ai = A2 = A, 
A 7^ XI, and therefore, the state space is not spanned by two linearly inde- 
pendent eigenvectors. 3 In this case, if u is an eigenvector of A, there is one 
solution with the usual form 



x\ (t) = e xt u. 

It can be verified that there is a second solution with the form 



X 2 (t) = e xt (tu + v). (2.24) 

The condition for (2.24) to be a solution is 

X 2 (t) = e xt [u + A (tu + u)] = Ax 2 {t) = e xt A{tu + v ) 
or, because e xt ^ 0, 

u = (A — XI) (tu + v). (2.25) 



3 The special case in which the (2 x 2) matrix A = A/, A real, sometimes called a bicritical node, 
has solutions that are half-lines from the origin, see figure 2.2(a). 
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Pre-multiplying by ( A — XI) we have 

(A- XI)u = (A- XI) 2 (tu + v). (2.26) 

But in view of the Cayley-Hamilton theorem mentioned above and given 
that, in this case, tr(H) = 2A and det(H) = A 2 , we have 

A 2 - 2XA + X 2 = (A - XI) 2 = 0. 

Hence, (2.25) is satisfied if u is an eigenvector of A and v solves the equation 
u = (A — X I)v. Then, by the superposition principle, the general solution 
can be written as 

x(t) = C\x\(t) + c 2 x 2 (i) (2.27) 

and 

X'(0) = C\U + c 2 u 

where ci and c 2 are two arbitrary constants depending on the initial condi- 
tions. Because v and u are linearly independent, they span the whole two- 
dimensional space. Thus, any arbitrary initial condition can be expressed 
as a linear combination of these two vectors. If the double eigenvalue A < 0, 
the solution x(t) — > 0 as t —>■ +oo because the lim ( _ +0O e xt t = 0, that is, 
the exponential contraction dominates the multiplicative effect of t. In the 
special case that the double eigenvalue A = 0, the solution takes the form 
x(t) = ( c\u + c 2 u) + c 2 fu, leading to an unstable equilibrium, divergence 
from which, however, grows linearly, not exponentially. 

Solutions have more complicated forms if multiplicities are greater than 
two. Broadly speaking, when A is a triple eigenvalue, the solution may 
include a term t 2 , when it is a quadruple eigenvalue, a term t 3 and so on. 
The general formula for repeated eigenvalues is the following: 

h n i - 1 

x(t) = E E (2.28) 

3 = 1 1=0 

where nj > 1 is the multiplicity of the jth eigenvalue (so tij = 2 means 
that A j = A J+ i); h < m is the number of distinct eigenvalues; J2 n j = m -> 
(j = 1 ,...,h). The vector- valued constants kji depend on the m arbi- 
trary initial conditions and the independent vectors. For further details see 
Blanchard et al. (1998), pp.282 ff. 

The transformation of A also takes a slightly different form, called a (real) 
Jordan canonical form, when there are multiplicities and the eigenvectors 
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do not span the entire space. The resulting matrix has diagonal blocks such 
as 

A 1 
0 A 

for a real, double eigenvalue A, or 

(Bj I 

V 0 Bj 

where Bj is as before, I is a (2 x 2) identity matrix and 0 is a (2 x 2) 
null matrix, for a double complex conjugate pair. For higher-order multi- 
plicities of eigenvalues the Jordan form may have more complicated diago- 
nal blocks. For a detailed discussion of the Jordan form see, for example, 
Hirsch and Srnale (1974). 



2.3 Continuous systems in the plane 

In applications, we very often encounter linear systems of two differential 
equations (or a differential equation of second order), either because the 
underlying model is linear or because it is linearised around an equilibrium 
point. Systems in two dimensions are particularly easy to discuss in full 
detail and give rise to a number of interesting basic dynamic configurations. 
Moreover, in practice, it is very difficult or impossible to determine exact 
values of the eigenvalues and eigenvectors for matrices of order greater than 
two. We therefore dedicate this section to continuous systems in the plane. 
The general form can be written as 



x 

y 




f ail 
\a2i 




(2.29) 



with x,y 6 R, a^- real constants. If det(vl) ^ 0, the unique equilibrium, for 
which x = y = 0, is x = y = 0. 

The characteristic equation is 

0 = det( Oll ~ A 1112 ,) 

V 021 CL22 — Ay 

= A 2 — (on + 022 ) A + (011022 — 012021 ) 

= A 2 — tr(vl) A + det(^4). 



and the eigenvalues are 



Ai, 2 = ^(tr(A) ± n/a) 



(2.30) 
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where A = ( [tr(^4)] 2 — 4det(.A)) is called the discriminant. For system 
(2.29) the discriminant is 

A = (on + 022) 2 — 4(an022 — ai2«2i) 

= (oil — 022 ) 2 + 4012021- 

The different types of dynamic behaviour of (2.29) can be described in 
terms of the two eigenvalues of the matrix A, which in planar systems 
can be completely characterised by the trace and determinant of A. In the 
following we consider nondegenerate equilibria (for which Ai and A 2 are both 
nonzero), unless otherwise stated. We distinguish behaviours according to 
the sign of the discriminant. 



CASE 1 A > 0 
form 



Eigenvalues and eigenvectors are real. Solutions have the 

x(t) = cie Xlt u^ + C 2 e X2t u ^ 
y(t) = c\e Xlt u^ + c 2 e X2t u' 2 ) - 



Off-equilibrium, the slope of the orbits in the (. x , y) plane asymptotically 
tend to a common, definite limit equal to 



t—>oc ax(t) 



= lim = lim 

£— >oo X t—> 00 



ciAie 






u{ ' + c 2 \2e X2t u 2 

c\\ie Xlt Ui^ 



(2) 



+ c 2 A 2 e A2 %9 1) 



(2.31) 



and either ci or C 2 is nonzero. Suppose that Ai > A 2 . Dividing throughout 
by e Xlt and taking the limit, we have 



r dy _ 

dx u w ' 
u i 



This is the slope of the eigenvector associated with the dominant eigen- 
value, that is, the largest eigenvalue, Ai. It defines the direction with which 
orbits asymptotically converge to the fixed point for t — > +00 or t — > — 00 . 

We have three basic subcases corresponding to figure 2.1(a), (b), (e), 
respectively (eigenvalues are plotted in the complex plane). 4 



(i) tr(^4) < 0, det(A) > 0. In this case, eigenvalues and eigenvectors 
are real and both eigenvalues are negative (say, 0 > Ai > A 2 ). The 
two-dimensional state space coincides with the stable eigenspace. The 

4 A special case is (let (.4) = 0 for which Ai = 0, A 2 = ir(.4). The solution corresponding to 
the zero eigenvalue is a constant {c\ e ^ 1 ' = ci) and the dynamics are determined by the sign 
of tr(A) as for cases (i) or (ii). In this case, however, there is no longer a unique equilibrium, 
but a continuum of equilibria whose locus is the straight line through the origin defined by 
y/x = —(an/ai2) = -(021/022)- 
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Fig. 2.1 Equilibrium types in the plane 



equilibrium is called a stable node, the term ‘node’ referring to the 
characteristic shape of the ensemble of orbits around the equilibrium. 

(ii) tr(yl) > 0, det(yl) > 0. In this case, eigenvalues and eigenvectors are 
real, both eigenvalues are positive (say, Ai > A 2 > 0) and the state 
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space coincides with the unstable eigenspace. The equilibrium is called 

an unstable node. 

(iii) det(.A) < 0. In this case, which implies A > 0 independently of the 
sign of the trace of A, one eigenvalue is positive, the other is nega- 
tive (say, Ai > 0 > A 2 ). There is, then, a one-dimensional stable and 
a one- dimensional unstable eigenspace and the equilibrium is known 
as a saddle point. All orbits starting off-equilibrium eventually di- 
verge from equilibrium except those originating in points on the stable 
eigenspace which converge to equilibrium. 

remark 2.4 Mathematically speaking, a saddle point is unstable. In eco- 
nomic literature, especially in rational expectations or optimal growth mod- 
els, we often find the concept of saddle point stability which, in the light 
of what we have just said, sounds like an oxymoron. The apparent con- 
tradiction is explained by considering that those models typically comprise 
a dynamical system characterised by a saddle point, plus some additional 
constraints on the dynamics of the system such that, given the initial condi- 
tion of some of the variables, the others must be chosen so as to position the 
system on the stable set. Therefore, the ensuing dynamics is convergence 
to equilibrium. 5 



When eigenvectors are real, as in the cases of nodes and saddle points, the 
eigenvalues and eigenvectors can be given interesting interpretations. To see 
this, let us assume that the first element of each eigenvector is nonzero and 
let us fix the arbitrary scale factor so that, for each vector, the first element 
is equal to one. Therefore we can write 

m = (l,uf) T 
u 2 = (1,4 2) ) T 



In the two-dimensional case, the expansion of (2.10) gives 



f an 

V«21 



au A 
«22 ) 




K 




whence 



1 (2) x 

a n + ai2o) = Xi 

, ( 2 ) \ ( 2 ) 

021 + 022 o) = A iU\ 



5 See for example, Magill (1977), p. 190; Brock and Malliaris (1989), p. 125; Turnovsky (1996), 
p. 137. 
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and the second element of the eigenvector is easily calculated as 



(2) A* — an 

u ) = 

0-12 



021 

Aj — 022 



i = 1,2. 



(2.32) 



Applying our discussion of section 2.2 to the simple two-dimensional case 

( 2 ) 

we can conclude that the quantities u\ ' are equal to the slopes of the straight 
lines emanating from the origin in the plane (x, y ) which are invariant under 
the law of motion defined by (2.29). From (2.29) the proportional rates of 
change of the two variables x, y are 



d\nx 

dt 

dlny 

dt 



x aux + auy , y 

- — — an + oi2 — 

XX X 

y a 2 ix + a 22 y x 

~ = =021 H 0 22 . 

y y y 



(2.33) 



Using (2.33) we can see that for initial conditions on one or the other of the 
lines defined by (2.32), we have 



x_y_ 

— — 'n 

x y 

that is, the proportional rates of change for both variables are constant 
and equal to the corresponding eigenvalue. The lines defined by (2.32) are 
sometimes denoted as balanced growth paths. On these paths the two 
variables evolve in time maintaining a constant ratio. These paths can be 
of special significance in applications. For example, in an economic context, 
if y denotes capital stock and x is output, along a balanced growth path we 
have a growing (or declining) output with a constant capital/output ratio. 
Alternatively, if y = x, and x is income, along the balanced growth path we 
have a constant proportional rate of growth (or decline) of income. 

Notice that on the balanced growth paths the value of the variables may 
increase or decrease according to the sign of the eigenvalues and the initial 
conditions. A positive proportional growth rate may result either from 
positive rate of change {x > 0) and positive level (x > 0) or from negative 
rate of change (x < 0) and negative level (x < 0). A negative proportional 
growth rate results from the rate of change and level being of opposite sign. 
Finally, from the discussion above, it follows that the balanced growth path 
corresponding to the largest eigenvalue (in algebraic terms) is attracting in 
the sense that almost all orbits (i.e. , all orbits except those starting on the 
other eigenvector) asymptotically converge to it in time. 
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CASE 2 A < 0 The eigenvalues and eigenvectors are complex conjugate 
pairs and we have 

(Ai, A2) = (A, A) = a ± 3/3 



with 

a = \ tr(-A) /3=^v /= A. 

As we have seen in the previous section, the solutions have the form 

x{t) = Ce at cos(/3f + <j>) 
y(t ) = Ce at sin(/3f + (f>) 



and the motion is oscillatory. If a 7^ 0 there is no strict periodicity in 
the sense that there exists no r such that x(t) = x(t + r). However, a 
conditional period can be defined as the length of time between two 
successive maxima of a variable, which is equal to 2n//3. The frequency 
is simply the number of oscillations per time unit, that is, (3/2n. The 
amplitude or size of the oscillations, depends on the initial conditions and 
e at (more on this point below). 

The complex conjugate eigenvectors can be calculated as follows. We once 
again use the degree of freedom and, assuming that the real (denoted Re) 
and imaginary (denoted Im) parts are nonzero, set = 1 + *, u ^ = 1 — %. 
Then, from (2.10), we obtain 



u 1 
02 



a — P — an a + /3 — on' T 

1 + l, h l 



012 



012 



A . a-/3-on .a + fi — anX 

1 -*, * 

V O12 012 / 



(2.34) 



From (2.34) we can construct two real and linearly independent vectors 



1 / . \ (a ^ ^ our 

e R = x( u i + U V = [ ) 

1 ( -a-/3 + an\T 

ei = x(o 1 - u 2 = ( - 1, 

2 v 0)2 > 

which span a plane that is invariant. Notice that in this particular case, the 
invariant plane coincides with the state space and with the stable or unstable 
eigenspace. This is no longer the case in higher-dimensional systems (see, 
for example, figure 2.4(b)). There are three subcases depending on the sign 
of tr(A) and therefore of Re(A) = a, see the corresponding figures in figure 
2.1(c), (d), (f), respectively: 
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(i) tr(A) <0, Re A = a < 0. The oscillations are dampened and the 
system converges to equilibrium. The equilibrium point is known as a 
focus or, sometimes, a vortex, due to the characteristic shape of the 
orbits around the equilibrium. In this case the focus or vortex is stable 
and the stable eigenspace coincides with the state space. 

(ii) tr(A) >0, Re A = a > 0. The amplitude of the oscillations gets larger 
with time and the system diverges from equilibrium. The unstable 
eigenspace coincides with the state space and the equilibrium point is 
called an unstable focus or vortex. 

(iii) tr(A) = 0, Re A = a = 0, det(A) > 0. In this special case we have a 
pair of purely imaginary eigenvalues. Orbits neither converge to, nor 
diverge from, the equilibrium point, but they oscillate regularly around 
it with a constant amplitude that depends only on initial conditions 
and a frequency equal to \/det(A)/27r. The centre eigenspace coincides 
with the state space and the equilibrium point is called a centre. We 
shall have more to say on this special case when we discuss nonlinear 
systems in chapters 3 and 5. 

CASE 3 A = 0 The eigenvalues are real and equal, Ai = A 2 = A. In 
this case, if A ^ XI, only one eigenvector can be determined, call it u = 
(u^, u^) T , defining a single straight line through the origin. Recalling the 
form of the second solution from (2.24), we can write the general solution 
as 

x(t) = (ciri*- 1 -* + C 2 'V^)e xt + tc 2 U^e xt 
y(t ) = ( c\u ^ + C 2 f ( ' 2 ' ) )e At + tC 2 U^e xt 

whence 

x(0) = C\U^ + C2'V^ 

y( 0) = ciu*- 2 -* + C 2 'C < ' 2 ^. 

The equilibrium type is again a node, sometimes called a Jordan node. An 
example of this type is provided in figure 2.2(b), where it is obvious that 
there is a single eigenvector. If we choose [x(0),y(0)] such that c\ ^ 0, 
C 2 = 0, we have, for every t > 0 

y(t) = y(0)_ = 

x(t) x(0) iii 1 ) 

thus the straight line defined by u is invariant. Along that line, the system 
asymptotically converges to the equilibrium (A < 0) or diverges from it 
(A > 0). However, the vector v does not define a second invariant straight 
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Fig. 2.2 Equilibrium types in the plane with repeated eigenvalue: (a) bicritical 

node; (b) Jordan node 

line. To see this, choose the initial conditions so that ci = 0, C 2 0. Then 
at t = 0 we have 

y(0) v 
x(0) ul 1 ) 

and at time t ^ 0 

y(t) = y{ Q) u {2) 

x(t) x(0) ut 1 ) 

which clearly is not time-invariant. 

Finally, if A = XI the equilibrium is still a node, sometimes called a 
bicritical node. However, all half-lines from the origin are solutions, giving 
a star shaped form (see figure 2.2(a)). 

The set of solution curves of a dynamical system, sketched in the state 
space, is known as the phase diagram or phase portrait of the system. 
Phase diagrams in two dimensions usually provide all the information we 
need on the orbit structure of the continuous-time system. The appendix 
explains the procedure to construct the phase diagram for typical cases of 
distinct eigenvalues and we refer the interested reader to it. 

Finally, figure 2.3 provides a very useful geometric representation in the 
(tr(H), det(H)) plane of the various cases discussed above. Quadrants III 
and IV of the plane correspond to saddle points, quadrant II to stable nodes 
and foci and quadrant I to unstable nodes and foci. The parabola divides 
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det (A) 




Fig. 2.3 Continuous-time dynamics in M 2 



quadrants I and II into nodes and foci (the former below the parabola, the 
latter above it). The positive part of the det(A)-axis corresponds to centres. 

The analysis for systems with n > 2 variables can be developed along 
the same lines although geometric insight will fail when the dimension of 
the state space is larger than three. In order to give the reader a broad 
idea of common situations we depict sample orbits of three-dimensional 
systems in figure 2.4, indicating the corresponding eigenvalues in the com- 
plex plane. The system in figure 2.4(a) has two positive real eigenvalues 
and one negative eigenvalue. The equilibrium point is an unstable saddle. 
All orbits eventually converge to the unstable eigenspace, in this case the 
plane associated with the positive real eigenvalues, and are captured by 
the expanding dynamics. The only exceptions are those orbits initiating 
on the stable eigenspace (defined by the eigenvector associated with the 
negative eigenvalue) which converge to the plane at the equilibrium point. 
Notice that the term saddle in M m refers to all the cases in which there 
exist some eigenvalues with positive and some with negative real parts. We 
use the term saddle node when eigenvalues are all real, saddle focus 
when some of the eigenvalues are complex. An example of the latter is 
represented in figure 2.4(b). The two real vectors associated with the com- 
plex conjugate pair of eigenvalues, with negative real part, span the stable 
eigenspace. Orbits approach asymptotically the unstable eigenspace defined 
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by the eigenvector associated with the positive eigenvalue, along which the 
dynamics is explosive. 



2.4 General solutions of linear systems in discrete time 

Let us consider now the linear, discrete-time system (2.5) 

x n+ i = Bx n x € M m . 

Once again we observe that x = 0 is an equilibrium solution and that for 
linear systems like (2.5) it is, generically, the only equilibrium. The general 
solution of (2.5) can be obtained following an argument entirely analogous 
to that used in section 2.2. If is a real, distinct eigenvalue of the (m x m) 
matrix B and Vi is the corresponding real eigenvector so that Bvi = KiVi, it 
can be verified that 



Xi(n) = k f Vi (2.35) 

is a solution of (2.5). At this point we ought to mention an important 
difference between solutions of continuous-time systems and solutions of 
discrete-time systems regarding the status of the invariant lines defined by 
the eigenvectors. The real eigenvectors associated with real eigenvalues for 
systems of difference equations are still invariant sets, that is, sequences of 
points initiating on these sets remain on the sets. However the invariant, 
continuous, straight lines are not solution orbits as they were in the con- 
tinuous case. Also recall from our discussion in chapter 1, section 1.2 the 
interesting peculiarity of discrete-time systems in the case of negative, real 
eigenvalues known as improper oscillations. In such a case, we may have 
the following interesting situation. Consider again system (2.5). Suppose 
we have one real eigenvalue k\ < 0, |«i| > 1 and that the initial condition 
is cc(0) = v\ (the eigenvector corresponding to «q). The solution is 

x{n) = k™v i (2.36) 

that is, the system stays forever on the straight line through the origin de- 
fined by v\ . At each iterate the system moves towards the equilibrium point 
x = 0, but always overshoots it, so that the distance from equilibrium be- 
comes ever larger in time, tending to oo as n —> +oo. Therefore, when there 
are negative real eigenvalues in discrete-time systems, great care should be 
applied in the use and, in particular, the interpretation of arrows in the 
phase diagram. 
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remark 2.5 Notice that if system (2.5) is characterised by improper oscil- 
lations then the system 

n ■ 1 — B~X n 

with B 2 = BB, will have the same stability properties, but no improper 
oscillations. In fact, if k is an eigenvalue of B, then p = k 2 is an eigenvalue 
of B 2 and \p\ = |ac| 2 . Then \p\ ^ 1 if |k| ^ 1. 

remark 2.6 If the (nonsingular) matrix B has an odd number of nega- 
tive eigenvalues (and therefore det(H) < 0), system (2.5) is said to be 
orientation-reversing (whereas for det(-B) > 0 it will be said to be 
orientation-preserving). The justification for these phrases can be read- 
ily seen in the simplest case of a one-dimensional system x n+ \ = G(x n ) = 
— x n , x n 6 R. The oriented interval [a, b] is mapped by G to an in- 
terval [G(a),G(6)j with reversed orientation. However, G 2 [a,b\ preserves 
orientation. 



G(b) G{a) 0 a b 



As usual, dealing with complex conjugate eigenvalues is more compli- 
cated. Suppose that we have a pair of eigenvalues of B 

(«i, • i ) = ('■'./• R i) = a .i ± ie 3 

with a corresponding pair of eigenvectors 

(vj,v j+ i) = = pj ± iqj 

where pj and qj are m-dimensional vectors. Then the pair of functions 

X M) = g + 

= r” [pj cos (ujn) - qj sin(wj-n)] 

x j+1 (n) = —(lijVj - /«"%) 

= r" \pj sin (ujjn) + qj cos(tUj?r)] , 

are solutions of (2.5). Here we have used the polar coordinate transforma- 
tions: 

a j = rj cos LOj 
9j = rj sin ujj 



(2.38) 
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and a well-known theorem due to De Moivre which states that 

(cosu; ± i sinu;) n = cos (can) ± isin(u;n). 

From (2.38) we have r = \J a 1 + 9‘ 2 . Then r is simply the modulus of the 
complex eigenvalues. The solutions in (2.37) can be simplified in a manner 
analogous to that used for continuous-time systems. If pp and q- 1 denote 
the fth elements of pj and qj , respectively, then in polar coordinates we have 

pf* = Cf* cos ((ftp) 
qf* = Cf sin 

where l = 1,2, Using again the trigonometric identities in (2.15), 

the solutions (2.37) can be re-written as m-dimensional vectors whose Zth 
elements have the form 

x j\n) = cos (ujjn + <f>^) 

Xj+i(n) = Cj l \jSm(ujjn + 4>^). 

Notice the similarities and differences with solutions to continuous-time sys- 
tems with complex eigenvalues in (2.16). 

Stable, unstable and centre eigenspaces are defined in a manner similar 
to that used for continuous-time systems, with the obvious difference that 
the discriminating factor for discrete-time systems is whether the moduli of 
the relevant eigenvalues are, respectively, smaller, greater or equal to one. 

Once again, assuming that we have m linearly independent solutions de- 
fined by (2.35) and (2.37), by the superposition principle the general solution 
of (2.5) can be written as a linear combination of the individual solutions 
as in (2.17), namely: 

x(n) = CiXi(n) + c 2 x 2 (n) H b c m x m (n) (2.39) 

where c* are constants depending on the initial conditions. 

As in the continuous case, a complication arises when eigenvalues are 
repeated. The general solution analogous to (2.28) is 

h n j~ l 

x ( n ) ~ kjin l tf 

j=i i=o 

where rij > 1 is the multiplicity of the j th eigenvalue, h < m is the number 
of distinct eigenvalues. 

Inspection of (2.35) and (2.37) indicates that if the modulus of any of the 
eigenvalues is greater than one, solutions tend to + or — oo as time tends 
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to plus infinity. On the contrary, if all eigenvalues have modulus smaller 
than one, solutions converge asymptotically to the equilibrium point. In 
the latter case, we say that the equilibrium is asymptotically stable, in the 
former case we say it is unstable. For distinct eigenvalues, the solution 
in the long run is determined by the term (or terms in the complex case) 
corresponding to the eigenvalue (or the pair of eigenvalues) with the largest 
modulus. When the dominant real eigenvalue or complex eigenvalue pair 
has modulus exactly equal to one, we have a weaker form of stability which 
we discuss in chapter 3. If that eigenvalue or eigenvalue pair is repeated, the 
equilibrium is unstable. As for continuous-time systems, there are criteria 
that can be applied to the coefficient matrix B to ensure stability. In this 
case, stability conditions guarantee that all eigenvalues have moduli inside 
the unit circle of the complex plane. There is also a compact form for the 
general solution of (2.5). We proceed as we did for the continuous-time case. 
Let X n be the matrix whose columns are the m solutions defined by (2.35) 
and (2.37). We can verify that 

X n = B n X o 

where X n = X(n) and Xq = X(0). Then 

X n — X n C — B Xq 

where c = (ci, C 2 , . . . , c m ) T are the constants of (2.39) and x n = x(n), 
xq = cc(0). To verify this consider that 

x n+ i = B n+l xo = BB n x 0 = Bx n 
that is, we have (2.5). 

As in the case of continuous-time systems, we often encounter discrete- 
time linear systems in two variables, either because the original system is 
linear or because we are studying a linearised system in two dimensions. 
We now turn to a detailed discussion of these. 



2.5 Discrete systems in the plane 

The discrete autonomous system analogous to continuous system (2.29) is: 



( Xn+l \ _ B f Xn \ _ / bn bl2 \ f Xn \ 

V Vn+ 1 ) V Vn ) V &21 b-22 ) V IJn ) 



(2.40) 



If the matrix (/ — B) is nonsingular, which we assume here, there exists a 
unique equilibrium point for (2.40), situated at the origin. The characteris- 
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tic equation is analogous to the continuous case, that is, 

k 2 — tr(£>) k + det(.B) = 0 

and the eigenvalues are 

« 1,2 = \ (tr (B) ± ^/[tr(B)] 2 -4det(5)j . 

The dynamics for the discrete system (2.40) are discussed case by case, 
again considering non degenerate equilibria (|«i|, |« 2 | 7^ 1): 

CASE 1 A > 0 The eigenvalues are real and solutions take the form 

x(n) = 

y{n) = + C 2 K 

(i) If |ki| < 1 and | « 2 ( < 1 the fixed point is a stable node. In discrete 
time this means that solutions are sequences of points approaching the 
equilibrium as n — > +oo. If k\, k? > 0 the approach is monotonic, 
otherwise, there are improper oscillations (see figures 2.5(a) and (c), 
respectively). In this case, the stable eigenspace coincides with the 
state space. 

(ii) If |ki| > 1 and |k 2 | > 1 the fixed point is an unstable node, that is, 
solutions are sequences of points approaching equilibrium as n — > — oo. 
If K \ , K 2 > 0 the approach is monotonic, otherwise, there are improper 
oscillations (as in figures 2.5(a) and (c), respectively, but arrows point 
in the opposite direction and the time order of points is reversed). In 
this case, the unstable eigenspace coincides with the state space. 

(iii) If | aci | > 1 and | AC 2 j < 1 the fixed point is a saddle point. No sequences of 
points approach the equilibrium for n ±oo except those originating 
from points on the eigenvectors associated with K\ or K 2 - Again, if 
«i, K 2 > 0 orbits move monotonically (see figure 2.5(b)), otherwise 
they oscillate improperly (see figure 2.5(d)). The stable and unstable 
eigenspaces are one-dimensional. 

CASE 2 A < 0 Then det(£>) > 0 and eigenvalues are a complex conjugate 
pair 

(ki, K 2 ) = (k, R) = a ± i9 

and solutions are sequences of points situated on spirals whose amplitude 
increases or decreases in time according to the factor r n where r = |cr±i#| = 




